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1.  Introduction 

In  this  paper,  we  report  simulations  of  turbulent  shear  flows  by 

numerical  solution  of  the  three-dimensional  Navier-Stokes  equations.  This 

) 

approach  provides  several  advantages  over  more  conventional  approaches. 

First,  the  complete  flow  field  is  obtained  at  all  times  so  that  detailed 

flow  characteristics  may  be  obtained  that  would  be  difficult  to  measure 

/ 

in  the  laboratory.  Second,  initial  conditions  can.  be  accurately  controlled 

t 

k 

so  that  their  effect  may  be  determined.  Third,  numerical  simulations  are 
convenient  experiments  to  assess  the  effect  of  various  physical  processes, 

.like  chemical  reactions,  on  turbulence  and  vice  versa.  Fourth,  the  simula- 
tions results  can  be  used  to  test  and  suggest  various  statistical  hypo- 
theses involved  in  turbn1  ci'ce  models  r.r.d  theories. 

The  present  work  is  an  extension  to  inhomogeneous  shear  flows  of  the 

three-dimensional,  homogeneous,  isotropic  turbulence  simulations  reported  by 

0rs2ag  & Patterson  (1972).  In  this  earlier  work,  simulations  using  up  to 
3 

(•32)  Fourier  modes  to  represent  each  velocity  component  were  performed 
at  microscale  Reynolds  number  20-50  (grid  Reynolds  numbers  5000-3000) , 
corresponding  to  moderate  Reynolds  number  grid  turbulence.  The  results  of 
these  homogeneous  turbulence  simulations  agreed  well  with  both  turbulence 

theory  (e.g.,  the  direct  interaction  approximation)  and  laboratory  experiments. 

• * 

Simulations  have  also  been  made  of  two-dimensional  'turbulence'  using  up  to 
(128)  Fourier  modes  (Herring  ct  al  1973)  with  successful  comparisons  with 
the  results  of  turbulence  theories.  These  successful  simulations  gave  impetus 
to  the  present  extension  to  shear  flows. 


In  tlie  homogeneous  turbulence  simulations,  it  was  observed'  t'hat  some 
important  features  of  the  flows  were  Reynolds  number  independent  (Herring  et 
al  1973),  so  that  there  is  basis  for  assuming  that  the  large-scale  features 
of  flows  simulated  at  huge  Reynolds  numbers.  It. is  also  reasonable  to  assume 
Reynolds-nunber- independence  of  large  scales  in  shear  flows  provided  initial 
and  boundary  conditions  are  Reynolds  number  independent.  It  seems  that 
Reynolds  number  dependencies  observed  in  laboratory  flows  may  be  ascribed  to 
variations  in  initial  or  boundary  conditions;  for  example,  Reynolds  number 
variations  in  turbulent  jets  seem  to  be  mainly  due  to  variations  in  inlet 
conditions  (like  boundary  layer  thickness)  with  Reynolds  number. 

* 

As  a first  application  of  our  shear  flow  turbulence  codes,  we  have 
simulated  the  monenturaless  wake  of  a self-propelled  body.  As  pointed  out 
by  Naudasc'ner  (1965),  the  momentumless  wake  bears  close  relationship  to  grid- 
generated homogeneous  turbulence,  since  both  are 'characterized  by  a very 

• • 

limited  region  in  which  there  is  significant  energy  transfer  from  mean  flow 
■to  turbulence  through  Reynolds  stresses. 

. , , ■ . . • ..  , , . i ‘ ‘ < 
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, In  Sec.. 2,  we  summarize  the  dynamical  equations  and  boundary  conditions 

employed  with  particular  emphasis  on  the  momentumless  wake  model.  In  Sec.  3, 
some  " ovel  aspects  of  the  numerical  approximation  are  discussed,  while  in 
Sec.  4,  some  results  are  presented  for  wake  turbulence.  Finally  in  Sec.  5,  we 
summarize  our  results  and  future  outlook. 

2 , Equation  of  Motion 

The  Navier-Stokes  equations  of  motion  for  an  imcompresslble  fluid  are  . 

» 

y (jc , t)xu  (jj  , t)-SlT  (jc , t)+vVZ v (jc , t) 

3t 


^•v(x,t)=0 
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where  y(x,t)  is  the  three-dimensional  velocity  field,  «5i(5S»t)=2  x y(?s,t) 

t » 

is  the  vorticity,  Tr(x,t)  » p + ^ is  the  pressure  he"'*,  p(x,t)  is  the  pressure, 
and  V is  the  kinematic  viscosity.  Eq.  (1)  is  written  in  rotation  form  to  fa- 
cilitate -numerical  solution  (Sec.  3). 


Boundary  conditions  require  more  discussion.  In  order  to  simulate  the 
momentumless  wake  of  a self-propelled  body  the  computational  domain  should 
include  a sizeable  region  of  potential  flow  both  upstream  and  downstream  of 
the  body  as  large  as  the  body  itself.  However,  this  would  be  very  wasteful  where 
most  of  the  computational  degrees  of  freedom  would  be  involved  in  resolving 
the  flow  outside  the  turbulent  wake.  Even  within  the  wake,  the  first  several  body 
diameters  downstream  are  of  adjustment  to  more  self-similar  conditions  down- 
stream. The  resolution  problem  is  a serious  one.  Presently  practicable  numeri- 
cal simulations  involve  at  most  order  of  10^  degrees  of  freedom  to  describe  the 
velocity  field,  there  would  remain  little  more  than  10^  degrees  of  freedom  to. 
do  the  flow  in  the  turbulent  wake.*  Clearly,  it  is  not  possible  to 
details  of  a turbulent  flow  with  so  little  resolution: 


In  order  to  avoid  the  resolution  problem  just  described,  we  simulate 
wake  flow  in  the  following  way.  We  isolate  a slab  in  a momentumless  turbulent 
wake  (Fig.  1)  and  follow  its  time  evolution  by  considering  it  enclosed  in 
a three-dimensional  box  as  shown  in  Fig.  2.  The  wake  axis  is  assumed  to 
be  along  the  xj  - axis  and  periodic  boundary  conditions  are  applied  at  Xj  = 

0,  and  y.^  ~ 0,  l^*  On  the  other  hand,  regid  free-slip  (no  stress)  or 
rigid  no-slip  boundary  conditions  are  applied  at  x^  - 0,  L^*  This  choice 
&f  boundary  conditions  in  the  x_-  direction  is  necessitated  by  the  face  that 
our  computer  codes  are  written  to  solve  also  the  Boussinesq  equation  of 
motion  of  a stratified  fluid  for  which  the  x^-ciirection  is  singled  as  the 
Urection  of  gravity  and  the  equations  of  shear  flow  in  a channel  with  rigid 


, l5,:' 
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top  and  bottom  boundaries. 


Initial  conditions  are  chosen  so  that  the  slab  of  turbulence  is  a 
realization  of  a section  of  a fully  developed  turbulent  wake,  i.e. 


v(x,0)  = v(x,0)  +v'(x,0)  (3) 

» Here  the  mean  defect  velocity  v(x,0)  is -chon.  .-f'to  be 

v(x,0)  = v *(r)  x (4) 

- d,  i 

^ 2 12  12 
where  Xj  is  the  unit  vector  in  the  x^-direction,  r^  =(X2_2~2^ ' +^X3~2L3^ 

and 

09  • 

/ v (r)r  dr=0 
0 d 

in  order  to  simulate  the  moraentumless  wake.  The  initial  fluctuating 

velocity  v* (x,0)  is  chosen  as  a realization  of  an  incompressible  random 
•» 

velocity  field  with  specified  local  energy  spectrum  and  turbulence.  Some 
'of  the  details  involved  in  the  construction  of  v*  are  explained  in  Sec.  3. 


! 

I 

i 


Numerical  solution  of  the  Navier-Stokes  equations  from  the  imposed 
initial  conditions  gives  the  time  evolution  of  the  simulated  cylindrical 
section  of  turbulence.  After  evolution  time  t,  the  results  are  interpreted 
as  a realization  of  a section  of  the  wake  flow  at  a distance  x^  = U^t  down- 
stream from  the  location  of  the  initial  wake  section,  where  Uq  is  the  body 
velocity.  In  other  words,  the  initial  flow  is  chosen  to  have  the  same  (or 
rather  similar)  statistical  properties  as  a section  of  a turbulent  wake  and 
the  time  evolution  of  the  flow  is  interpreted  as  the  downstream  variation  of 
the  wake.  The  model  we  solve  numerically  is  statistically  homogeneous  along 
| the  wake  axis  x.  but  nonstationary  in  time;  the  wake  in  the  frame  of  a uni- 
formly  moving  body  is  statistically  stationary  in  time  but  is  inhomogeneous 
i Xj.  It  is  asserted  that  the  Galilean  transformation  Xj  =Ugt  relates  the 
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numerical  and  physical  experiments. 


3.  Numerical  Methods 

3.1  Free-slip  boundaries 

If  the  boundary  conditions  at  x^  = 0,  are  no-stress  (free-slip)  while  • 


those  at  x^  **  0, 

and  x2  c Qv  L2  are  periodic  then 

t 

3v 

3V2  ' 

*=  « V3  on  x3  = 0,L3;  v(x1-)mL1,x2+nL2,x3)  = v(x) 

(6) 

In  this  case,  the  velocity  field  can  be  expanded  in  the  Fourier  series 


va(2*t>  " 


y 


y 

x— r ua  (k , t ) exp  [ 2tt  i ( k^/Lj+k^  /L2  ^ ^ 


cos  7ik0x-/L-  Ct  = 1,2 
3 3 3» 

sin  nk^x^/L^,  a - 3 


(7) 


where  the  indicated  summations  are  over  integer  k^,  k2,  k^.  For  the  .simulations 

reported  below,  the  cutoffs  are  chosen  as  = K2  = 16,  = 32,  so  that  the 

3 

spectral  representations  (7)  each  involve  about  (32)  independent  degrees  of 
freedom  (Fourier  amplitudes). 

We  have  used  the  expansions  (7)  in  two  kinds  of  numerical  approximation 
to  the  Navier-Stokes  equations.  In  the  special  method  (Orszng  1971a),  equations 
for  u(k)  are  derived  by  substituting  (7)  into  (1),  multiplying  the  result  by 
exp[-2iri(k^x^/LjL  + k2x2/L2)]cos [nk^x^/L^]  for  a = 1,  2 (and  by  the  same 
expression  with  the  cosine  replaced 'by  sine  if  a = 3),  and  finally  integrating 
the  result  over  the  box  0 £ xft  < LQ.  The  resulting  equations  are,  after  elirainat 
of  the  pressure  by  means  of  the  incompressibility  constraint  (2) , 


•-.TCil’K  
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. • 3ua(k,t)  x~_  _ _ . _2  'J*  _ 

~af  + vk  “a^5  = •lk6(Wv/k  > uB(E>,)uY(3’t) 

“Kp  < Pp»  9p  < Kp 

where  k » 2nk  /L  for  a = 1,2  and  k0  = ttk./L.,  and 
a a a 3 3 3. 


ua(k,t)  = 


1 (k^ , k2 , | k^ | * t )^  a ~ 1,2 
2i  ^3u3^i*^2*  ,tXa“  ^ 


(8) 


(9) 


Numerical  solution  of  (8)  is;  accomplished  using  fast  Fourier  transform 
methods  to  evaluate  the  convolution  sums,  leapfrog  time  differencing  on  the 
nonlinear  terms,  and  Crank-Nicols'on  (implicit)  time  differencing  on  the 

3 

viscous  terms.  Overall,  18  Fourier  transforms  on  (32)  points  must  be  per- 
formed each  time  step  (Orszag  1971a). 

In  the  pseudospectral  method  (Fox  and  Orszag  1973) , the  expansions  (7) 
are  used’ as  an  interpolatory  tool  to  evaluate  the  derivatives  appearing  in 
(1).  "Grid1'  points  x^  — o^j^/2K^,  — u,  ...»  zK  — x for  a — 1,2,  x^  — 


era  a’  Ja 


4 - 0 } # * j K ~tc  introduced  ^nd  $2ri2s  C^/  srs  vsssd  to  ovElust? 


derivatives  as,  for  example, 
J3 


to  (x,t)  = E E E [ik  u9(k,t)-ik,u. (k, 

3-  IkJ^  |k2|<K2  0<k3<K3  21- 


O) 


(10) 


• exp[ik^x^+ik2x2 3 cos  k^Xg 


The  final  result  of  this  procedure  is  to  give  a set  of  equations  for  the 
spectral  amplitudes  u(k),  defined  by  (9),  that  are  identical  to  (8)  except  for 
the  replacement  of  the  convolution  sums  by  similarly  truncated  sums  with 
pQ  + qa  = ka  + 2K  or  p a + qa  = ka  (Orszag  1971a).  The  additional  terms 
entering  the  sums  in  (8)  are  usually  called  "aliasing"  terms.  It  has  been 
shown  (Fox  and  Orszag  1973)  that  the  differences  in  the  results  obtained  by 
the  present  pseudospectral  method  and  :he  spectral  method  are  generally 
•>gligible  so  long  as  either  method  gives  an  accurate  solution  of  the  equations 
otion.  Since  the  pseudospectral  method  may ‘be  implemented  in  only  9 


t 


!\  Fourier  transforms  over  (32) ^ points  per  time  step,  it  is  roughly  a factor 

2 more  efficient  than  the  spectral  method  and  so  has  been  used  for  most  of 
\ the  simulations  reported  below. 

| ' Both  the  spectral  and  pseudospectral  methods  have  been  programmed  for 

| a CDC  7600.  The  programs  involve  double  buffering  of  data  between  small 

r 

i core  memory,  large  core  memory,  and  disks.  The  spectral  method  requires 

! 

j about  6‘  s per  time  step  on  the  CDC  7600,  while  the  pseudospectral  code 
& • 

t 

v requires  only  3.2  s per  time  step,  both  for  the  cutoffs  K,  **  K9  = 16,  K,  *»  32. 

$ 

Many  of  the  critical  internal  loops  of  the  program  are  written  in  assembly 
5 language  to  avoid  inefficiencies  attributable  to  the  Fortran  compiler, 

5 

I although  further  improvements  in  the  code  should  permit  a speedup  of  nearly 


Besides  the  speed  advantage  of  the  pseudospectral  method  over  the  spectral 
method,  the  former  has  the  great  advantage  Hint-  it  nnnl-ips  with  but  the  most 
minor  of  modifications  to  problems  involving  more  complicated  physics,  like 
chemical  reactions  or  radiation.  Since  the  expansions  (7)  are  used  merely  as 
> an  interpolator^  tool  in  the  evaluation  of  derivatives,  they  may  be  similarly 
used  in  the  evaluation  of  these  more  complicated  effects.  Nevertheless,  the 
pseudospectral  method  shares  all  the  advantages  of  the  spectral  method  with 
regard  to  accuracy  and,  especially,  efficiency  improvement  over  finite- 
difference  schemes  (Orszag'fc  Israeli  1974).  The  expression  of  (1)  in  rotation 
form  is  useful  since  it  gives  pointwise  energy  conservation. 

3.2  No-slip  boundaries 

| With  no-slip  boundary  conditions  applied  nt  x^  = 0,  L^»  the  velocity  satisfies 

* * 

1 v *»  0 on.  x^  B 0,Lg>  vCxjdml^^+nl^.x^)  = v(x)  (11) 

instead  of  (6).  It  is  no  longer  appropriate  to  use  the  Fourier  expansions 
'7),  not  just  because  v^  a V2  ° 0 at  * 0,  1.^,  but  rather  because  imposition 
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of  a Fourier  series  representation  of  the  dependence  would  induce  Gibbs 
phenomena  at  the  boundaries  and  result  in  slow  convergence  of  the  Fourier 
series  (Orszag  1971a). 

No-slip  boundaries  are  best  treated  by  Chebyshev  expansion  in  , 

(Orszag  1971b).  The  velocity  field  is  represented  as 

u(k,t)  exp[ik^x^+ik2X2 
'3 

■ (12) 

where  the  nth  degree  Chebyshev  polynomial  T^(x)  is  defined  by  Tn(cosO)  = 
cos  n0.  It  may  be  shown  that,  if  v(x)  is  smooth,  the  series  (12)  or  any  of 
its  termwise  derivatives  do  not  exhibit  Gibbs  phenomena  at  the  boundaries. 
Equivalently,  the  series  (12)  converges  faster  than  algebraically  as  K -*». 

Notice  that  the  boundary  conditions  (11)  must  be  imposed  as  constraints 
on  (12). 

One  advantage  of  (12)  is  that  it  leads  to  pseudut>pecLi.dl  dppj.uAiu»atiuuS 
to  (1),  (2)  that  are  very  similar  in  form  to  that  following  from  the  Fourier 
series  (7).  In  particular,  the  pseudospectral  equations  with  (12)  may  be 
implemented  in  9 Fourier  transforms  per  time  step,  and  are  but  slightly  less 
efficient  than  for  free-slip  boundaries.  With  cutoffs  ^ - 16  and  = 32, 

our  code  requires  4 s per  time  step.  In  this  latter  code,  tine  differencing 

is  done  by  Adams-Bashforth  differencing  (Lilly  1965)  on  the  nonlinear  terms 

% 

(to  avoid  instability  that  would  result  from  use  of  leapfrog  because  of  the 
stability  induced  by  the  boundary  conditions)  and  Cxank-Licolson  on  the  viscous 
terms.  The  pressure  computation  is  done  by  using  (2)  to  get  an  equation 
tridiagonal  in  the  Chebyshev  index  k,  for  the  pressure  and  diagonal  in 
and  k2«  Solution  of  the  resulting  tridiagonal  system  accounts  for  most  of 
the  additional  time  required  by  the  rigid  boundary  code  over  the  free-slip 


iTfc  ([2x3-L3l  /L3) 


v(x,i)  = 


E E 2 


l<Ki 


ik2l<K2 


,0<k3< 


ndnrv  mAo. 


* 


Methods  based  on  the  spectral  expansion  (12)  have  an  important  advantage 

over  difference  methods,  in  addition  to  the  advantages  they  enjoy  when 

periodic  of  free-slip  boundary  conditions  are  applied.  The  "grid”  points 

for  the  pseudospectral  method  based  on  (12)  are  x^  - 1 (1  + cos 

for  = 0,  K^»  so  that  the  effective  resolution  near  the  walls  x^  = 0,1^ 

is  Ax^  ^ L^/K^ • *n  fact,  if  there  is  a boundary  layer  of  thickness  6 along 

either  x^  = 0 or  x^  = L^j  it  may  be  shown  (Orszag  & Israeli  1974)  that  it 

1/2 

■ is  sufficient  to  take  'u  3(L^/S)  to  achieve  better  than  17.  accuracy 
in  the  boundary  layer.  In  effect,  the  Chebyshev  polynomial  expansion  gives 
a highly  nonuniform  grid  near  the  boundaries.  This  behavior  is  particularly 
appropriate  for  the  study  of  channel  flows,  etc.  where  thin  boundary  layer 
! * are  apt  to  develop. 

I 

3.3  Initial  conditions 

•In  the  wake  model  introduced  in  Sec.  2,  initial  conditions  of  the  form 
of  a "cylinder"  of  turbulence  are  required.  The  mean  defect  velocity  (4)  may 
be  imposed  arbitrarily,  but  the  fluctuating  component  v’ (x)  must  satisfy  the 
incompressibility  constraint.  This  is  done  by  writing 


v*  (x)  « V x A*  (x)  (1 

where  the  vector  potential  A* (x)  is  chosen  to  achieve  the  desired  local 

energy  spectrum  and  turbulence  intensity.  It  suffices  to  choose  A* (x)  to  be 

> ^ — 


i of  the  form 


U 

A’(x)  = ri(r)r  B(X)  (14) 

^ •• 

•where  the  turbulence  intensity  function  I(r)  is  a nonrandom  function  of 
only  the  distance  from  the  wake  axis  x^  and  the  fluctuation  component  B(x) 

| is  chosen  as  a realization  of  a homogeneous,  isotropic  random  field  with 

specified  isotropic  energy  spectrum,  as  done  by  Orszag  & Patterson  (1972).  If 


(W*srt 
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the  intensity  function  is  chosen  to  vanish  outside  a radius  from  the  wake 
axis,  the  resulting  turbulent  velocity  field  is,  by  (13),  nonturbulent 

s' 

outside  the  cylinder  of  radius  r^  centered  on  the  wake  axis. 

In  summary,  our  technique  cf  imposing  the  initial  conditions  allow 
arbitrary  mean  velocity  profile,  turbulence  intensity  profile,  and  local 

i 

turbulence  energy  spectrum. 

4.  Momentumless  Wake 

The  wake  model  of  Sec.  2 and  the  numerical  methods  of  Sec.  3 permit 
simulation  of  the  turbulent  wake  of  a self-propelled  body.  Free-slip  boundary 
conditions  are  applied  at  x^  = 0,  and  the  spectral  cutoffs  = K2  = 
e 32  are  used.  The  following  choice  of  initial  parameters  was  made: 

Li  « L2  *=  L3  = 2,  vd(r)  = Vq  sin  (r/r1)/(r/r1)  where  r1  - 1/4,  I(r)  « 
oax[l-r  /4.5»0],  and  turbulence  energy  spectrum  E(k)  - Ak  u.n.p(-Bk  ) [cf.  Gvszag 
& Patterson  1972].  The  initial  conditions  arc  chosen  to  match  as  closely 
as  possible  the  results  of  Naudascher  (1965)  and  Wang’  (1965)  at  a location 
4 body  diameters  behind  a self-propelled  disk  in  a wind  tunnel.  In  the 
numerical  simulations,  the  viscosity  is  chosen  to  be  V = .005,  so  that  the 
Reynolds  number  of  the  simulation  is  roughly  13,000  (run  2a),  and  V = .003 

Vith  Reynolds  number  roughly  22,000  (run  3).  In  comparison,  the  laboratory 

% 

experiments  of  Naudascher  and  Wang  were  run  at  Reynolds  numbers  of  roughly 
55,000. 

Some  measure  of  the  degree  of  complication  of  the  flow  that  is  simulated 
*s  given  by  Figs.  3&4.  In  Fig.  3,  contours  of  the  run  3 axial  mean  velocity 
arc  plotted  at  t = .3,  corresponding  to  about  7 diameter.'  downstream  from 
the  body  upon  making  identification  of  body  velocity  by  the  ratio  max[vd(r)] /IP 
4t  the  station  at  x/D  - 4 with  the  experimental  results.  The  axial  mean  is 
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determined  by  averaging  over  the  x^-direction.  In  Fig.  4,  contours  of  the 
run  3 x^-component  of  vorticity  are  plotted  at  t = .8. 

In  Fig.  5,  we  compare  the  axial  variation  of  maximum  turbulent  intensity 
with  the  results  of  Wang  and  Naudascher.  Here  Uq  and  t^  are  determined  as 
explained  above  by  correspondence  with  the  data  at  the  station  at  x = 4D, 
while  Xq  - 2D  according  to  Naudascher  and  Wang.  It  is  apparent  from  these 
results  that  the  present  simulations  are  in  substantial  agreement  with  the 
laboratory  results  of  Naudascher  and  Wang,  at  least  over  the  limited  downstream 
range  of  the  present  experiments. 

In  Fig.  6,  we  compare  the  radial  variation  of  axial  mean-square  turbulent 
intensity  in  the  laboratory  experiments  and  the  numerical  calculations ' The 
curve  labelled  t - 0 shows  the  initial  distribution,  while  that  labelled  t = .5! 
shows  the  resulting  distribution  for  run  2a  at  about  6D  downstream  from  the 
body.  Again,  the  agreement  with  the  experimental  results  is  satisfactory. 

It  is  apparent  from  the  simulation  results  of  this  section  that  numerical 
1 ’ simulation  of  turbulent  shear  flows  is  well  within  present  computational 

i . • ■ ■ • 

; capabilities.  However,  the  process  of  extracting  useful  information  about 

' shear  flows  from  numerical  simulations  is  still  in  its  infancy.  The  most 

significant  problem  with  the  present  simulations  is  the  scarcity  of  points 

\ 

in  the  axial  direction -with  which  to  compute  averages.  With  just  32  points 

in  the  longitudinal  direction,  statistical  errors  are  more  large  and  either 

% 

multi-time  step  information  or  ensembles  must  be  employed  to  improve  the 
statistics.  This  situation  should  be  contrasted  with  the  case  for  homogeneous, 
j isotropic  turbulence  where  space  averages  suffice  to  give  statistical  results 
generally  to  within  5%  (Orszag  & Patterson  1972). 


7 * *s*'7'*iK*f***~**+*+~~~ 


5.  Conclusions 

The  present  simulations  of  turbulent  shear  flows  are  but  a first  step 
toward  proper  understanding  of  the  basic  mechanisms  and  dynamics  of  these 

flows.  Detailed  comparisons  and  tests  are  presently  being  made  between 

>f  ! 

various  turbulence  modelling  hypotheses,  laboratory  experiments,  and  the 

i 

present  simulations.  As  time  and  the  art  of  numerical  simulation  progress, 
simulations  like  the  present  ones  should  be  expected  to  fulfill  more  and 
more  the  need  of  a laboratory  workhorse.  Out  present  simulations  provide 
ample  evidence  for  this.  The  computer  codes  we  have  written  and  sufficiently 
general  to  study  channel  flows  like  plane  Poiseuille  and  plane  Couette  flows, 
as  well  as  momentumless  and  momentumfull  wakes,  both  in  homogeneous  and 
stratified  fluids.  Experimental  set-up  of  this  variety  of  shear. flows,  not- 
withstanding specifying  the  form  of  the  shear  and  turbulence  profiles  that 

may  be  imposed,  would  be  a herculean  task.  On  the  other  hand,  the  computer 

« 

simulations  can  handle  all  these  cases. 

In  future  work  on  the  momentumless  wake,  we  shall  report  on  longer  simulation 
runs  now  underway  and  on  techniques  to  improve  the  statistics  of  the  results. 

In  order  to  improve  the  choice  of  initial  conditions,  simulations  runs  are  made 
in  which  the  pseudorandom  initial  conditions  imposed  as  described  in  Sec.  3.3 
are  let  to  evolve  for  about  10  body  diameters  downstream  and  then  are  reapplied, 

amplified  in  excitation,  at  a virtual  upstream  point  in  order  to  begin  the 

* 

calculation  anew.  This  approach  avoids  the  difficulty  of  modelling  imprecisely 
known  laboratory  experiments.  With. the  initial  conditions  imposed  as  in  Sec.  3.3, 
the  wake  is  very  sensitive  to  the  initial  turbulence  level  relative  to  the  mean 

t 

shear,  while  much  of  this  dependence  is  avoided  by  the  present  technique. 
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Fig.  2 Spatial  box  enclosing  the  section  of  wake  (0,  L^)  is  in 
Fig.  1..  Periodic  boundary  conditions  are  applied  at  the 
sidewalls  and  either  free-slip  or  no-slip  conditions  are 
applied  at  the  top  and  bottom. 


ng  (Naudascher)  1965 


Axial  variation  of  maximum  turbulent  intensity 
Cross  circles :• results  of  Wang  and  Naudascher. 
Squares : Run  2 Circles:  $un  3* 


